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Abstract 

Variational  assimilation  is  used  to  combine  velocity  and  sea-surface  height  anomaly  (SSHA)  measurements  with  a 
system  of  dynamics  to  estimate  the  seasonal  flow  through  the  Tsushima/Korean  Strait  for  the  summer,  autumn,  and 
winter  seasons  of  1999-2000.  The  velocity  measurements  are  from  two  lines  of  moored  acoustic  Doppler  current 
profilers  (ADCPs)  spanning  the  Tsushima/Korean  Strait  just  north  and  south  of  Tsushima  Island  and  the  SSHA 
measurements  are  from  the  TOPEX  altimeter.  The  dynamics  are  the  linear,  time-independent,  shallow-water  equations 
and  are  forced  by  winds  from  the  Navy  Global  Ocean  and  Atmospheric  Prediction  System.  A  weighted  least-squares 
technique  is  used  to  determine  the  seasonal  flow  fields  that  simultaneously  minimize  the  weighted  residuals  of  the  two 
data  sets  and  the  system  of  dynamics.  The  weights  are  based  on  expected  errors,  allowing  the  assimilation  system  to  put 
more  emphasis  on  the  components  of  data  and  dynamics  that  are  known  more  accurately. 

Earlier  studies  show  that  the  flow  through  the  Tsushima  Strait  is  barotropic  throughout  most  of  the  year.  The  region 
just  offshore  of  the  southeast  coast  of  South  Korea,  however,  has  been  identified  as  being  subject  to  strong  baroclinic 
processes  during  autumn.  In  an  attempt  to  verify  this  observation,  best  estimates  of  the  seasonal  flow  fields  are 
computed  by  assimilating  the  two  data  sets  with  both  a  single-layered  barotropic  system  of  equations  and  a  2.5-layered 
baroclinic  system  of  equations.  The  minimized  residuals  of  the  barotropic  system  of  equations  reveal  that  the  solutions 
satisfy  these  equations  within  expected  errors  throughout  most  of  the  domain.  Just  offshore  of  the  southeast  coast  of 
South  Korea,  however,  the  barotropic  momentum  equation  residuals  exceed  3  times  the  expected  error  during  autumn, 
therefore  implying  that  the  barotropic  approximation  is  not  valid  at  this  location.  The  2.5-layered  baroclinic 
momentum  equation  residuals  within  this  region  are  smaller  than  the  barotropic  residuals,  suggesting  the  existence  of 
baroclinic  processes. 
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1.  Introduction 

The  Tsushima/Korea  Strait  (Fig.  1)  is  the 
primary  inlet  of  seawater  and  advected  heat  into 
the  Japan/East  Sea.  The  majority  of  water  that 
flows  into  the  Tsushima  Strait  originates  from  the 
Kuroshio  and  its  properties  are  distinguished  by 
being  relatively  warm  and  saline.  Water  enters  the 


strait  through  two  different  locations.  The  major¬ 
ity  of  incoming  water  enters  between  Kyushu  and 
Cheju  Island  and  is  called  the  Tsushima  Current 
(TSC).  The  remaining  inflow  comes  in  between 
Cheju  Island  and  South  Korea  and  is  called  the 
Cheju  Warm  Current  (CWC).  Tsushima  Island 
divides  the  incoming  water  into  the  eastern  and 
western  channels  (Fig.  1).  The  number  of  distinct 
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Fig.  1.  Map  of  the  Tsushima  Strait  region  where  contours  indicate  bathymetry  (m)  and  the  black  dots  represent  the  discretized  grid  at 
which  best  estimates  of  vertically  integrated,  seasonally  averaged  velocity  and  seasonally  averaged  SSH  are  computed.  Also  displayed 
on  this  map  are  the  positions  of  the  bottom-mounted  ADCPs  (gray  stars). 
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currents  flowing  out  of  the  Tsushima/Korean 
Strait  and  into  the  Japan/East  Sea  is  observed  to 
vary  from  one  to  three.  The  exiting  current  that  is 
most  persistently  observed  is  the  Near  Shore 
Branch  (NSB).  This  current  is  an  extension  of 
the  TSC,  originates  from  the  eastern  channel,  and 
flows  northward  along  the  coast  of  Japan  (Kim 
and  Legeckis,  1986).  The  water  exiting  the  western 
channel  typically  forms  the  East  Korean  Warm 
Current  (EKWC),  which  flows  northward  along 
the  Korean  continental  slope.  Observations  have 
been  made,  however,  that  show  that  the  EKWC 
does  not  exist  during  winter  and  spring  (Kim  and 
Legeckis,  1986;  Cho  and  Kim,  2000).  A  third 
current  is  occasionally  observed  at  the  exit  of  the 
Tsushima/Korean  Strait  and  is  referred  to  as  the 
Off  Shore  Branch  (OSB)  (Kim  and  Legeckis, 
1986).  This  current  forms  as  the  EKWC  begins 
to  turn  northward  and  flow  between  the  coast  of 
Korea  and  the  Ullung  Basin.  A  branching  occurs 
on  the  downstream  side  of  the  strait,  and  a  mass  of 
water  flows  eastward  along  the  southern  edge  of 
the  Ullung  basin.  Continuing  alongside  the  edge  of 
this  basin,  the  mass  of  water  turns  northeastward, 
forms  the  OSB,  and  enters  into  the  Japan/East  Sea 
flowing  parallel  to  the  NSB. 

The  origin,  branching,  and  properties  of  the 
flow  through  the  strait  have  been  observed  to  have 
substantial  seasonal  and  interannual  variations. 
The  average  transport  through  the  strait  is  roughly 
3.5  Sv  (1  Sv  =  106  m3/s)  (Miita  and  Ogawa,  1984). 
This  transport,  however,  has  been  observed  to 
vary  from  5.6  Sv  in  September  to  1.2  Sv  at  other 
times  (Isobe  et  al.,  1994).  In  addition,  since  this 
strait  is  only  about  250  km  across  with  depths 
typically  less  than  150  m,  the  flow  through  the 
strait  is  subject  to  a  strong  constraining  influence 
by  the  coastlines.  The  Rossby  Radius  of  deforma¬ 
tion  is  approximately  18  km  (Cho  and  Kim,  1998). 
Strong  seasonal  variations  in  surface  forcing,  such 
as  wind  stress  and  heat  flux,  also  play  an 
important  role  in  causing  the  environment  within 
the  strait  to  change.  Since  the  bathymetry  within 
the  strait  is  shallow,  surface  forcing  is  able  to 
penetrate  and  influence  a  majority  of  the  water 
column.  During  winter,  when  the  wind  stress  is  at 
its  maximum  and  the  transport  is  at  its  minimum, 
the  water  within  the  strait  is  well-mixed  and  the 


flow  is  observed  to  be  barotropic.  Whereas,  during 
summer  and  autumn,  when  surface  forcing  is 
relatively  weak  and  the  transport  swift,  a  strong 
thermocline  at  30-40  m  depth  is  observed  (Jacobs 
et  al.,  2001b).  Isobe  (1995)  and  Cho  and  Kim 
(1998)  hypothesize  that  baroclinic  processes  are 
prevalent  during  these  warmer  seasons,  especially 
at  the  exit  of  the  western  channel. 

In  an  attempt  to  verify  these  baroclinic  pro¬ 
cesses  and  to  estimate  the  seasonal  flow  through 
the  strait,  an  approach  is  used  to  obtain  a  solution 
that  best  fits  prescribed  dynamics  and  measure¬ 
ments  in  a  weighted  least  squares  sense  (varia¬ 
tional  assimilation).  Two  sets  of  measurements 
contribute  to  the  assimilation.  The  first  set  includes 
velocity  measurements  from  an  array  of  moored 
acoustic  Doppler  current  profilers  (ADCPs)  de¬ 
ployed  during  summer,  autumn,  and  winter  of 
1999-2000.  The  second  set  is  seasonal  sea-surface 
height  (SSH)  change  inferred  from  TOPEX 
altimeter  data.  A  description  of  these  two  sources 
of  data,  including  the  processing  and  preparation 
for  assimilation,  is  given  in  the  next  section.  The 
dynamics  that  the  solution  is  constrained  to  follow 
is  described  in  Section  3.  Two  sets  of  numerical 
dynamics  are  used  to  test  the  accuracy  of  the 
barotropic  approximation  in  the  strait.  These 
dynamics  are  the  linear,  time-independent,  shal¬ 
low-water  equations  and  boundary  and  smoothing 
constraints  for  both  a  single-layered  barotropic 
system  and  a  2.5-layered  baroclinic  system. 

A  solution  is  constructed  that  simultaneously 
minimizes  the  dynamic  and  data  residuals,  where 
the  residuals  are  the  misfits  of  the  numerical 
representation  of  the  dynamic  and  data  caused 
by  a  particular  solution.  These  residuals  are 
combined  into  a  cost  function,  such  that,  given 
an  estimate  of  the  flow  field,  the  cost  function 
provides  a  single  scalar  value  that  represents  the 
total  misfit  to  dynamics  and  data.  The  optimal 
solution  is  determined  by  performing  a  direct 
variational  minimization  of  this  cost  function, 
which  amounts  to  solving  a  weighted  least  squares 
problem  (Le  Dimet  and  Talagrand,  1986).  In  order 
to  properly  construct  a  solution,  the  weighting  of 
the  data  and  dynamical  equations  must  be 
correctly  specified  in  the  cost  function.  This  is 
done  by  estimating  the  level  of  weight  (inverse  of 
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the  expected  error  covariance)  for  each  compo¬ 
nent.  Certain  methods  are  available  to  determine 
expected  error  covariances.  The  expected  errors  of 
the  measurements  can  be  determined  directly, 
whereas  the  error  levels  of  the  dynamical  equa¬ 
tions  are  more  difficult  to  quantify.  Therefore,  it  is 
necessary  to  test  a  range  of  expected  error  levels  to 
ensure  that  the  solution  is  not  overly  sensitive  to 
the  specified  weights.  This  assimilation  method, 
including  the  specification  of  the  expected  errors,  is 
described  in  Section  4. 

The  seasonal  circulation  fields  resulting  from  the 
assimilation  of  both  ADCP  and  TOPEX  data  with 
the  barotropic  system  of  equations  are  presented 
and  discussed  in  Section  5.  To  demonstrate  the 
impact  that  these  two  data  types  have  on  the 
solution  and  its  cost  function,  results  from 
assimilating  three  different  configurations  of  just 
ADCP  data  are  also  presented  and  discussed. 
Since  these  solutions  are  least-squares  fits  to 
observations  and  dynamics,  they  satisfy  neither 
exactly.  Examination  of  the  dynamic  residuals 
reveals  where  and  during  which  seasons  the 
solution  deviates  significantly  from  the  prescribed 
equations  of  motion.  These  results  suggest  that  the 
flow  through  the  Tsushima/Korean  Strait  is 
mostly  barotropic.  The  barotropic  momentum 
equation  residuals,  however,  are  significantly  high 
during  autumn  near  the  exit  of  the  western 
channel.  Since  these  residuals  far  exceed  expected 
values,  it  is  implied  that  baroclinic  processes  might 
be  an  important  contributor  in  this  region. 
Support  for  this  hypothesis  is  provided  by 
comparing  these  residuals  with  the  residuals 
resulting  from  using  a  system  of  2.5-layered 
baroclinic  equations.  Analysis  of  this  comparison 
reveals  that  the  solutions  satisfy  the  baroclinic 
momentum  equations  better  than  the  barotropic 
momentum  equations,  especially  in  the  region  and 
during  the  season  where  baroclinic  processes  are 
suspected. 

2.  Assimilated  data 

In  this  study,  two  types  of  data  are  used  in  the 
estimation  of  seasonal  circulation  fields  (Fig.  2). 
The  first  data  set  is  comprised  of  seasonally 


averaged  vertically  integrated  velocities  derived 
from  moored  current  meter  data  (vertically  in¬ 
tegrated  velocity  is  hereinafter  referred  as  trans¬ 
port).  The  second  set  of  data  is  seasonal  SSH 
differences  inferred  from  TOPEX  range  measure¬ 
ments.  Even  though  these  two  data  sets  represent 
different  features  and  have  different  levels  of 
accuracy,  by  using  variational  assimilation  the 
relative  accuracy  is  accounted  for  in  the  weighting 
used  in  the  least-squares  process  (discussed  in 
Section  4.3). 

2.1.  Velocity  measurements  from  ADCP  moorings 

The  velocity  measurements  are  obtained  from 
13  bottom  mounted  ADCP  moorings  deployed  by 
the  Naval  Research  Laboratory.  Due  to  the  high 
concentration  of  fishing  and  trawling  that  occurs 
throughout  the  strait,  it  has  been  difficult  in  the 
past  to  deploy  moorings  over  long  periods  of  time. 
In  order  to  overcome  this  problem,  trawl  resistant 
bottom  mounts  (called  Barnys  after  their  barnacle¬ 
like  shape)  were  used.  The  exposed  side  of  these 
dome-shaped  mounts  is  relatively  smooth  in  order 
to  minimize  snagging  by  fishing  nets  and  lines,  and 
the  instrumentation  is  enclosed  safely  within  the 
dome.  Each  mooring  was  equipped  with  a  RD 
Instruments  Workhorse  ADCP  operating  at 
300  KHz  (Perkins  et  al.,  2000).  These  moorings 
were  deployed  in  May  of  1999  and  spanned  the 
width  of  the  Tsushima  Strait  just  northeast  and 
southwest  of  Tsushima  Island  (Fig.  1).  An  addi¬ 
tional  mooring  was  deployed  in  October  of  1999  in 
the  western  channel  (mooring  Cl). 

The  moorings  were  recovered  in  March  of  2000, 
and  all  data  was  available  and  used  except  from 
moorings  N1  and  Cl  during  summer  and  autumn, 
and  mooring  S5  during  winter.  The  collected  data 
consists  of  both  u  (eastward)  and  v  (northward) 
components  of  velocity  with  individual  samples 
containing  an  RMS  error  of  1  cm/s.  The  temporal 
resolution  of  these  samples  is  15  min  at  stations  S5, 
N2,  and  N3,  and  30 min  at  the  other  stations.  The 
vertical  resolution  is  2m  at  stations  SI,  S2,  and 
N2,  and  4  m  at  the  other  stations.  With  bottom- 
mounted  ADCPs,  however,  measurements  close  to 
the  surface  are  not  useable  due  to  contamination 
of  the  acoustic  signal  by  surface  interference.  The 
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Fig.  2.  Sea-surface  height  anomaly  (SSHA)  measurements  from  the  TOPEX  altimeter  (colored  contours)  and  vertically  integrated 
velocity  measurements  from  the  bottom  mounted  ADCPs  (red  vectors)  are  seasonally  averaged  for  (A)  summer,  (B)  autumn,  and  (C) 
winter  of  1999-2000.  The  SSHA  is  an  interpolation  of  TOPEX  measurements  to  all  discretized  grid  points  (Fig.  1)  within  50  km  of 
TOPEX  groundtracks  (numbered  gray  lines).  The  interpolation  used  is  a  linearly  weighted  average  based  on  distance  to  the  data 
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distance  from  the  surface  at  which  this  interference 
occurs  varies  with  instrument  deployment  depth. 
The  relatively  deep  ADCPs  (about  140  m)  were 
able  to  measure  currents  to  within  20  m  of  the 
surface,  whereas  the  ADCPs  in  the  shallower 
regions  (about  100  m)  were  able  to  measure 
currents  to  within  10  m  of  the  surface  (Jacobs  et 
al,  2001a). 

Tidal  variability  within  the  strait  accounts  for 
20-60%  of  the  velocity  variability  depending  on 
the  depth  and  mooring  location  (Teague  et  al., 

2001) .  Since  this  study  is  concerned  with  just  the 
nontidal  variability,  tidal  currents  and  currents 
near  the  tidal  semidiurnal  and  diurnal  frequencies 
are  removed  from  the  ADCP  data  using  a  low-pass 
filter  with  a  40-h  cutoff  frequency  (Teague  et  al., 

2002) .  The  measurements  of  u  and  v  at  each 
location  are  vertically  integrated  over  either  the 
entire  water  column  for  the  barotropic  system  or 
each  of  two  layers  for  the  baroclinic  system.  These 
transports  are  then  rotated  43°  in  the  counter¬ 
clockwise  direction  so  that  the  u  component  is  in 
the  along-strait  direction  to  match  the  discretized 
grid  (Fig.  1).  Finally,  these  measurements  are 
seasonally  averaged  over  summer,  autumn,  and 
winter  of  1999-2000  (Fig.  2).  The  seasonal  periods 
are  chosen  by  examining  the  total  transport 
through  the  strait.  Distinct  transitions  in  the 
transport  are  observed  between  seasons.  The 
transition  between  summer  and  autumn  occurs 
rapidly,  whereas  the  transition  between  the  au¬ 
tumn  and  winter  transport  is  relatively  slower  and 
requires  about  a  month  to  occur  (Teague  et  al., 
2002).  The  dates  that  are  best  separated  by  these 
transitions  are:  summer  (5/15/1999-8/15/1999), 
autumn  (8/15/1999-11/15/1999)  and  winter  (12/ 
15/1999-3/15/2000). 

2.2.  Sea  surface  height  measurements  from  the 
TOP  EX  altimeter 

Seasonal  differences  in  SSH  are  inferred  from 
processed  TOPEX  data  that  are  in  terms  of  sea- 
surface  height  anomaly  (SSHA).  Since  the  TOPEX 
orbit  and  range  measurements  are  known  accu¬ 
rately,  3.5  and  3.2  cm,  respectively,  the  sea  level 
relative  to  the  Earth  reference  ellipsoid  can  be 
determined  with  an  accuracy  of  4.7  cm  (Fu  et  al., 


1994).  Unfortunately,  this  measured  sea  level  gives 
little  useful  information  in  regards  to  the  SSH 
associated  with  dynamic  motion.  In  order  to 
obtain  this  information,  the  geoid  height  must  be 
subtracted  from  the  sea  level.  The  geoid  height  is 
the  distance  from  the  reference  ellipsoid  to  the 
geopotential  surface  (the  ocean  surface  if  it  were  at 
rest  with  no  forcing).  However,  since  this  geoid 
height  is  not  known  accurately,  a  time-average  of 
all  sea-level  measurements  made  from  the  first  8 
years  of  the  TOPEX/Poseidon  mission  is  used 
instead.  This  mean  sea  level  contains  both  the 
geoid  and  the  mean  height  due  to  ocean  circula¬ 
tion.  Therefore,  deviations  from  this  mean  are 
SSH  anomalies  (SSHA).  These  processed  TOPEX 
data  are  a  product  of  the  Naval  Research 
Laboratory  and  are  described  in  Jacobs  et  al. 
(2002). 

There  are  six  TOPEX  groundtracks  within  the 
vicinity  of  the  Tsushima/Korean  Strait  (Fig.  2). 
The  TOPEX  data  that  are  used  include  all 
available  SSHA  measurements  during  the  same 
time  period  that  the  moorings  were  deployed 
(cycles  245-276,  where  each  cycle  encompasses 
approximately  10  days).  For  each  cycle,  a  weighted 
average  of  all  SSHA  measurements  at  groundtrack 
points  within  50  km  of  each  model  grid  point  is 
used  to  compute  the  SSHA  on  the  model  grid.  The 
weight  decreases  linearly  from  one  at  0  km  from 
the  groundtrack  point  to  zero  at  50  km  distance. 
The  gridded  SSHA  data,  which  includes  976  grid 
points,  is  seasonally  averaged  over  the  same  time 
periods  that  are  used  with  the  velocity  measure¬ 
ments  (Fig.  2).  Directly  assimilating  this  gridded 
SSHA  data,  however,  can  result  potentially  in 
significant  error.  To  avoid  this  error,  seasonal 
differences  in  SSHA  are  computed  between 
autumn  and  summer,  and  autumn  and  winter, 
and  used  in  the  assimilation.  By  using  seasonal 
SSHA  differences,  the  reference  level  needed  to 
infer  SSH  becomes  irrelevant.  There  are  two 
drawbacks  to  this  approach;  the  first  is  that  all 
three  seasonal  solutions  must  be  solved  simulta¬ 
neously,  thus  increasing  the  computational  bur¬ 
den.  The  second  drawback  is  that  there  are  no 
direct  SSH  measurements  within  the  assimilation, 
therefore  allowing  the  SSH  solutions  to  be  biased 
by  a  value  that  is  constant  in  space  time. 
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Table  1 


Description  of  the  variables  and  constants  found  in  the  barotropic  shallow  water  equations  (Eqs.  (l)-(3)) 


Element 

Name 

Units 

Value 

U 

Vertically  integrated,  seasonally  averaged  velocity  in  the  along-strait  direction 

m2/s 

Solution  variable 

V 

Vertically  integrated,  seasonally  averaged  velocity  in  the  across-strait  direction 

m2/s 

Solution  variable 

n 

Seasonally  averaged  sea  surface  height 

m 

Solution  variable 

f 

Coriolis  parameter 

l/s 

Prescribed  function 

9 

Gravity 

m/s2 

9.81 

H 

Depth 

m 

Prescribed  function 

Horizontal  eddy  viscosity  coefficient 

m2/s 

1000 

Tx 

Wind  stress  in  the  along-strait  direction 

N/m2 

Prescribed  function 

A 

Wind  stress  in  the  across-strait  direction 

N/m2 

Prescribed  function 

Po 

Average  density 

kg/m3 

1024.24 

CD 

Bottom  drag  coefficient 

0.002 

I«bI 

Average  magnitude  of  bottom  velocity 

m/s 

Prescribed  function 

3.  Dynamics 


3.1.  Single-layered  barotropic  equations  of  motion 


The  equations  of  motion  used  to  describe  the 
flow  within  the  single-layered  barotropic  system 
are  based  on  the  shallow-water  equations.  The 
governing  physics  in  this  system  are  the  conserva¬ 
tion  of  mass  and  momentum.  These  equations  are 
a  simplification  of  the  fundamental  Eulerian 
equations  of  motion  that  emphasize  motions  with 
large  horizontal  length  scales  compared  to  the 
vertical  scale  (Lx,Ly$>  Lz).  In  addition  to  the 
hydrostatic  assumption,  Boussinesq  and  linear 
approximations  are  also  applied.  These  simplified 
equations  of  motion  are  Reynolds  averaged  over 
seasonal  time  scales  so  that  the  time  rate  of  change 
of  all  variables  is  neglected  and  the  equations 
represent  the  dynamics  of  the  seasonally  averaged 
flow  fields, 


dU 

dx 


dV  n 
+  -r—  —  0, 
dy 


jy 


CD|fiB|C/ 

H 


a  2u 


a 2u  d2v 


dx2  9 y1  3x3^ 


Po 


(2) 


-fU-gHp-  +  KH 
dy 


a2  v 


d2U 


dx2  ^  dy2  ^  dxdy 


CD\uB\V 

H 


,_Z 

Pi)' 


(3) 


where  U,  V  and  fj  are  the  transport  components 
and  SSH,  respectively. 

Each  of  the  variables  and  constants  in  Eqs. 
(l)-(3)  is  listed  along  with  their  units  and  values  in 
Table  1.  The  Coriolis  parameter  (/)  is  considered  a 
function  of  latitude.  The  horizontal  eddy-viscosity 
coefficient  (Kf)  is  highly  dependent  on  the 
problem  and  can  have  a  wide  range  of  values. 
Therefore,  the  value  of  Kn  is  estimated  by 
considering  the  values  and  grid  resolution  used  in 
many  existing  models.  The  average  density  (p0)  is 
estimated  using  the  Modular  Ocean  Data  Assim¬ 
ilation  System  (MODAS),  which  is  described  in 
Section  3.2.  The  value  chosen  for  the  bottom  drag 
coefficient  (CD)  is  a  commonly  used  value  for  the 
sea  floor  (Csanady,  1984).  The  average  magnitude 
of  bottom  velocity  (|wB|)  is  estimated  for  each 
season  individually  by  averaging  all  ADCP  velo¬ 
city  measurements  within  the  bottom  10  m  of  the 
water  column.  The  bathymetry  (H),  which  is  a 
product  of  the  Naval  Research  Laboratory,  has  a  1/ 
12°  resolution  and  is  interpolated  to  the  grid  using  a 
bicubic  spline  fit.  The  wind  stresses  (tv.  and  xy)  come 
from  the  Navy  Global  Ocean  and  Atmospheric 
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Prediction  System  which  has  a  spatial  resolution  of 
1°  and  a  temporal  resolutions  of  6h.  This  data  set 
is  interpolated  to  the  model  grid  points  using  a 
bicubic  spline  and  seasonally  averaged  over  the 
same  dates  as  the  measurements. 


3.2.  2.5-Layered  baroclinic  equations  of  motion 


The  baroclinic  system  of  equations  differs  from 
the  barotropic  system  in  that  density  is  no  longer 
considered  constant,  therefore  allowing  the  density 
to  provide  internal  pressure  gradients.  For  exam¬ 
ple,  by  allowing  the  density  to  be  a  function  of  x, 
y,  layer,  and  season,  the  pressure  gradient  term  in 
Eq.  (2),  becomes 


gH 


drj 

0x 


gHl(hz*LsJi+m„°A%i 

\  2pj  ox  ox )  2pl  ox 


(4) 


for  the  first  layer  of  the  baroclinic  system,  where  hx 
is  the  reference  thickness  of  layer  1,  is  the  SSH, 
and  fj2  is  the  deviation  of  the  interface  between 
layers  1  and  2. 

The  baroclinic  system  of  equations  is  split  into 
three  layers  with  the  bottom  layer  being  at  rest 
(If  =  Vi  =  0).  The  only  influence  that  the  bottom 
layer  has  on  the  equations  of  motion  is  that  the 
deviation  of  the  interface  between  it  and  layer  2 
(>f3)  affects  the  thickness  and  the  limits  of  vertical 
integration  of  the  second  layer.  Therefore,  the 
equations  of  motion  are  only  applied  to  the  first 
two  layers,  resulting  in  seven  unknown  variables 
(U i,  V\,  fju  U2,  Vj,  rj2,  and  fj3)  that  need  to  be 
solved  for  at  each  grid  point.  The  reference 
thickness  of  layer  1  (/q)  is  chosen  such  that  it 
approximately  contains  the  mixed  layer  that  is 
apparent  during  summer  and  autumn.  Through 
examination  of  the  MODAS  data,  the  average 
depth  of  this  mixed  layer  is  roughly  h\  =  45  m.  The 
reference  thickness  of  layer  2  (h 2  =  130  m)  is 
chosen  such  that  the  only  places  where  the  bottom 
layer  exists  within  the  gridded  domain  are  the  deep 
basin  that  is  west  of  Kyushu  and  the  Ullung  basin 
(Fig.  1). 

The  densities  and  density  gradients  in  the 
baroclinic  momentum  equations  are  computed 
independently  of  the  inverse  problem  using 
MODAS.  MODAS  is  one  of  the  present  US  Navy 


standard  tools  for  production  of  static,  bimonthly, 
three-dimensional  grids  of  temperature  and  sali¬ 
nity  using  all  available  historical  observations. 
Further  details  regarding  this  system  are  provided 
in  Fox  et  al.  (2002).  The  grid  resolution  of  the  data 
set  produced  by  this  tool  for  the  Tsushima  Strait 
region  is  1/8°  along  the  horizontal  plane  and  every 
10  m  up  to  50  m  depth  and  then  every  50  m  deeper 
than  50  m.  Densities  are  computed  from  this 
climatological  data  set  using  the  international 
equation  of  state  of  seawater,  1980  (Pond  and 
Pickard,  1983).  The  density  data  at  each  depth 
interval  are  interpolated  to  the  grid  points  using  a 
bicubic  spline  fit.  Finally,  by  seasonally  and 
vertically  averaging  the  interpolated  data  set  for 
each  layer  and  season,  the  functions  px  and  p2  are 
computed.  The  horizontal  gradients  of  px  and  p2 
are  computed  using  first-order  accurate  finite 
differencing. 


3.3.  Physical  constraints 


In  addition  to  continuity  and  momentum,  two 
additional  constraints  are  included  in  both  systems 
of  dynamics.  The  first  constraint  forces  the  cross¬ 
shore  component  of  flow  to  be  zero  at  the  solution 
grid  points  that  border  land.  This  constraint  is 
enforced  through  satisfaction  of  the  following 
equation: 


Uhx  +  Vhy  —  0, 


(5) 


where  hx  and  hy  are  the  normal  components  in  the 
along  and  cross-strait  directions,  respectively,  of 
the  vector  perpendicular  to  the  coastline.  These 
normal  vector  components  are  calculated  using  the 
normalized  gradients  of  bathymetry. 


0///0X 

dH/dx)2  +  (QH/dy)2  ’ 

i  =  dH/Qy 

^(dH/dx)2 +  (dH/dy)2’ 


(6) 


For  the  case  of  the  2.5-layered  system,  this 
boundary  constraint  is  applied  to  both  layers.  All 
of  the  domain  boundary  grid  points  that  border 
water  are  considered  open  and  have  no  constraint. 
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The  second  constraint  is  smoothing  and  its 
purpose  is  to  remove  features  that  can  not  be 
resolved  by  the  discretized  grid  or  the  measure¬ 
ments.  In  order  to  prevent  these  unresolved  small 
features  in  the  solution,  smoothing  is  imposed  by 
constraining  the  Laplacian  plus  a  second-order 
cross  derivative  of  each  of  the  solution  variables  to 
equal  zero. 


d2U  &U  d2U 

(7) 

0x2  '  8y2  '  0x0y 

d2V  d2V  d2V 

(8) 

dx2  8y2  dxdy 

S2>7  _0 

0x2  dy2  0x0y 

(9) 

Without  this  constraint,  the  solution  could  be 
filled  with  many  sub-observation-scale  features  not 
resolved  by  the  measurement  arrays.  For  example, 
given  a  fine  numerical  grid  with  a  coarse  set  of 
measurements,  a  small  feature  that  satisfies  the 
dynamics  can  be  placed  between  the  measurement 
positions  without  affecting  the  solution  at  the 
measurement  locations.  The  distance  between  the 
ADCP  moorings  is  roughly  25  km,  and  previously 
observed  currents  have  scales  of  28-37  km  (Katoh, 
1993).  Therefore,  the  ADCP  measurements  can  be 
expected  to  resolve  the  dominant  features  within 
the  vicinity  of  the  moorings.  For  the  case  of  the 
2.5-layered  system,  the  smoothing  constraints  in 
Eqs.  (7)  and  (8)  are  applied  to  the  upper  two 
layers,  and  Eq.  (9)  is  applied  to  all  three  layers. 

4.  Method 

4.1.  Domain 

Solutions  to  the  seasonally  averaged  transports 
and  SSH  are  computed  on  a  discretized  grid  (Fig. 
1).  Since  the  Tsushima/Korean  Strait  is  not 
aligned  in  either  the  meridional  or  zonal  direction, 
the  domain  is  rotated  43°  in  the  counterclockwise 
direction  so  that  the  local  x-axis  of  the  domain  is 
in  the  along-strait  direction.  This  is  done  to 
minimize  the  number  of  grid  points  required  to 
represent  the  region  adequately  while  maintaining 


relative  ease  in  the  computer  implementation  of 
the  numerical  solution.  The  domain  consists  of  62 
grid  points  in  the  along-strait  direction,  and  33  in 
the  cross-strait  direction  with  a  grid  spacing  of 
Ax  =  Aj>=10km.  The  resulting  arrangement 
of  grid  points  is  somewhat  irregular  due  to 
landmasses  and  islands.  Grid  points  are  removed 
from  the  surface  layer  if  they  are  on  land  or  in 
water  with  depths  less  than  10  m.  As  for  the  second 
and  third  layer,  grid  points  are  removed  if  the 
depth  is  less  than  the  thickness  of  the  layer  above 
it.  Also,  grid  points  that  are  not  a  part  of  at  least 
four  contiguous  water  grid  points  in  both  the 
along  and  cross-strait  directions  are  neglected 
(four  consecutive  grid  points  are  required  in 
order  to  perform  a  second-order  accurate  finite 
difference  operation  at  a  boundary  point). 
Application  of  this  masking  to  the  domain 
results  in  1545,  1450,  and  179  grid  point  locations 
in  the  surface,  middle,  and  bottom  layers,  respec¬ 
tively,  at  which  solutions  are  sought  for  each  of  the 
three  seasons.  Therefore,  the  total  number  of 
unknowns  that  need  to  be  solved  for  in  the 
barotropic  and  baroclinic  systems  are  13,905  and 
27,492,  respectively. 

4.2.  Weighted  least  squares 

The  solution  is  defined  as  the  set  of  seasonally 
averaged  transports  and  SSH  that  minimizes  the 
cost  function,  which  contains  the  weighted  resi¬ 
duals  of  the  discretized  equations  of  motion, 
physical  constraints,  and  measurements.  The 
equations  of  motion  (Eqs.  (l)-(3))  and  smoothing 
constraints  (Eqs.  (7)-(9))  are  discretized  at  each 
grid  point  for  each  of  the  three  seasons.  The 
derivative  terms  in  these  equations  are  discretized 
using  second-order  accurate  finite  differencing. 
The  boundary  constraint  (Eq.  (5))  is  used  only  at 
the  grid  points  that  border  land.  Contributions  to 
the  cost  function  due  to  residuals  between  the 
ADCP  measurements  and  solution  are  applied  at 
the  grid  points  closest  to  each  of  the  ADCP 
mooring  locations;  whereas,  the  SSHA  difference 
measurements  are  applied  at  the  locations  of  the 
interpolated  SSHA  (Fig.  2).  Since  the  SSHA 
difference  data  involve  more  than  one  season,  the 
equations  for  all  three  seasons  must  be  solved 
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simultaneously,  therefore  resulting  in  roughly 
30,000  equations  for  the  barotropic  system. 

The  barotropic  system  of  independent  equations 
(30,000)  provides  a  number  of  equations  greater 
than  the  number  of  unknowns  (13,905).  The 
solution  to  this  problem  is  determined  using  the 
weighted  least  squares.  To  begin,  the  system  of 
linear  discretized  equations  can  be  written  in 
matrix  form 

Ax  =  b,  (10) 

where  A  is  the  coefficient  matrix  which  contains 
the  discretization  of  the  left-hand  side  of  the 
momentum,  continuity,  boundary,  smoothing,  and 
measurement  equations.  The  state  vector  (x) 
contains  the  13,905  unknown  variables  for  which 
a  solution  is  sought  (U,  V,  and  fj  at  all  grid  points), 
and  the  forcing  vector  (b)  consists  of  the  right- 
hand  side  of  all  the  equations.  For  example,  the  x- 
momentum  equation  (Eq.  (2))  applied  at  a 
particular  grid  point  is  contained  in  one  row  of 
A,  and  the  wind  stress  at  this  grid  point  specifies 
one  element  of  b. 

All  of  these  discretized  equations,  however,  have 
at  least  some  error  associated  with  them. 

Ax  -  b  =  e  (11) 

The  task  at  hand  is  to  determine  the  solution 
(x)  that  minimizes  the  residuals  (e).  To  accomplish 
this,  the  residuals  are  combined  into  a  total 
cost  function  in  a  weighted  least-squares 
fashion, 

Cost  Function  =  J(x)  =  (Ax  —  b)TW(Ax  -  b), 

(12) 

where  W  is  a  weighting  matrix  that  places 
emphasis  on  the  equations  that  have  smaller 
expected  errors.  The  minimum  of  this  cost  func¬ 
tion  is  determined  by  setting  the  first  variation 
with  respect  to  the  state  equal  to  zero. 

J(x  +  <5x)  -  J(x)  =  0.  (13) 

By  taking  the  limit  of  this  equation  as  <5x  ->•  0, 
the  weighted  least-squares  formulation  is  pro¬ 
duced. 

ATWAx  =  ATWb.  (14) 


For  this  study,  an  iterative  conjugate  gradient 
technique  is  used  to  solve  Eq.  (14).  This  conjugate 
gradient  method  is  a  relaxation  method  that  begins 
with  an  initial  state  of  all  zeros  and  seeks  out  the 
final  state  by  moving  in  search  directions  that  are 
constructed  by  the  conjugation  of  the  residuals 
(Ax  -  b).  The  conjugate  gradient  method  is 
beneficial  since  the  search  directions  do  not  have 
to  be  stored  in  memory;  once  a  new  search 
direction  is  computed,  the  old  one  is  forgotten. 
In  addition,  the  conjugate  gradient  method  is  able 
to  take  advantage  of  the  sparseness  of  the  normal 
matrix  (ATWA)  in  Eq.  (14).  The  normal  matrix, 
however,  is  ill-conditioned  since  the  ratio  between 
the  largest  and  smallest  eigenvalue  is  very  large, 
and  at  least  one  eigenmode  is  known  to  be  ill- 
posed.  This  ill-posed  eigenmode  is  present  due  to 
the  lack  of  direct  SSH  measurements,  thus  causing 
the  overall  magnitude  of  the  SSH  solution  to  be 
dependent  on  the  initial  state.  Even  though  this 
inverse  problem  is  ill-conditioned  and  contains  an 
ill-posed  eigenmode,  best  estimates  of  the  flow  field 
can  be  computed  and  provide  reliable  information 
(Wunsch,  1996;  Thacker,  1992). 

To  test  the  consistency  of  the  solution,  dozens  of 
experiments  are  performed.  In  each  of  these 
experiments,  a  different  initial  state  is  constructed 
using  a  random  number  generator  and  then  used 
to  initialize  the  conjugate  gradient.  As  expected, 
the  magnitude  of  the  SSH  solution  varied  sub¬ 
stantially  from  experiment  to  experiment. 
However,  the  spatial  and  temporal  gradients  of 
the  SSH  solution  and  the  transport  solution 
never  varied  by  more  than  15%  of  the  original 
solution.  The  accuracy  of  the  final  solution  is 
dependent  on  the  number  of  iterations  through 
the  conjugate  gradient.  Theoretically,  for  a 
problem  that  is  well-conditioned,  this  routine 
should  be  able  to  converge  to  a  solution  in  a 
number  of  iterations  that  is  equivalent  to  one  to 
three  times  the  length  of  the  state.  However,  since 
the  problem  posed  in  this  study  is  ill-conditioned, 
which  means  there  are  many  solutions  that  are 
close  to  the  best  solution,  the  conjugate  gradient 
requires  many  times  this  number  of  iterations  to 
navigate  to  the  absolute  minimum  cost  within  the 
machine’s  precision  error.  This  is  discussed  further 
in  Section  5.2. 


S.R.  Smith  et  al  /  Deep-Sea  Research  II 52  (2005)  1763-1783 


1773 


4.3.  Weighting  matrix 

One  of  the  most  important  aspects  of  any 
inversion  or  statistical  estimation  scheme  is  the 
covariance  errors  or  the  weight  used  in  the  cost 
function.  Once  the  idealized  physical  processes  are 
embedded  within  the  dynamical  equations,  addi¬ 
tional  constraints  applied,  and  measurements 
taken,  some  thought  and  care  must  be  placed  into 
deriving  the  appropriate  expected  error  (<er>). 
These  errors  should  take  into  account  the 
neglected  terms  in  the  dynamic  equations,  errors 
in  forcing,  errors  in  the  measurements,  and 
representativeness  errors  (Daley,  1993).  After  the 
solution  is  obtained,  the  residuals  in  the  dynamical 
equations  are  compared  to  the  estimated  expected 
errors  to  determine  if  the  solution  is  in  agreement 
with  the  initial  error  estimates.  If  the  residuals  of 
the  dynamics  are  greater  than  the  a  priori  error 
estimates,  then  there  must  be  additional  contribu¬ 
tors  to  the  residuals  that  were  not  originally 
considered. 

The  weighting  matrix  that  is  used  in  this 
assimilation  is  a  diagonal  matrix  whose  elements 
consist  of  the  reciprocal  of  the  expected  variance 
l/<er2>  for  each  corresponding  dynamical  equa¬ 
tion.  This  means  that  the  cross-covariance  of  the 
expected  errors  are  neglected.  For  simplicity,  it  is 
also  assumed  that  the  expected  variance  for  each 
of  the  equations  of  motion,  constraints,  and 
measurements  are  the  same  at  all  grid  points  and 
during  all  seasons. 

In  order  to  maintain  the  correct  proportion  of 
weights  within  the  total  cost  function  so  that  they 
are  independent  of  the  number  of  grid  points  or 
measurements,  the  weight  for  each  dynamical 
equation  is  divided  by  its  number  of  discretized 
equations  for  each  season.  For  example,  the 
continuity  weight  is  defined  as 

^Continuity  —  l/C^Continuityl^T"  )  Continuity)*  (15) 

where  ^continuity  is  the  number  of  points  at  which 
the  continuity  equation  is  applied,  n  is  also  equal 
to  the  number  of  grid  points  for  the  momentum 
equations  and  smoothing  constraints.  For  coastal 
boundary  constraints  and  measurements,  n  is  the 
number  of  grid  points  that  border  land  and  the 
number  of  measurements  that  are  assimilated  for 


each  season,  respectively.  This  definition  of  the 
weights  assures  that  each  contribution  to  the  cost 
function  can  be  modeled  as  a  chi-squared  variable 
with  an  expected  value  of  one  and  is  of  comparable 
magnitude.  For  example,  if  the  contributions  to 
the  cost  function  are  not  divided  by  n ,  a  doubling 
of  resolution  would  place  greater  weight  by  a 
factor  of  four  on  the  dynamical  equations. 

Optimally,  the  expected  error  associated  with 
each  of  these  dynamical  equations  should  include 
all  of  the  error  resulting  from  neglected  terms, 
approximations,  and  discretization.  Since  it  is 
nearly  impossible  to  account  for  all  of  the  sources 
of  error  within  each  dynamical  equation,  the 
estimation  of  the  expected  errors  entails  some 
error.  Therefore,  the  expected  error  associated 
with  each  set  of  dynamical  equations  is  approxi¬ 
mated  by  only  considering  the  primary  source  of 
error. 

The  primary  source  of  error  in  the  continuity 
equation  (Eq.  (1))  is  a  result  of  neglecting  the  time 
rate  of  change  of  SSH  (dfj/dt).  This  term  is 
estimated  by  seasonally  averaging  the  TOPEX 
data  at  each  of  the  data  points  within  the  region. 
Then,  at  each  data  point  the  SSHA  differences 
between  consecutive  seasons  are  computed. 
Averaging  the  magnitudes  of  these  differences 
results  in  a  mean  seasonal  SSHA  difference  of 
A fj  =  17.67  cm.  The  estimated  value  for  the  ne¬ 
glected  term  dfj/dt,  and  hence  the  estimated 
expected  error  of  continuity  is  given  in  Table  2 
along  with  the  expected  errors  associated  with  the 
other  sets  of  dynamics  and  data. 

In  the  derivation  of  the  momentum  equations 
(Eqs.  (2)  and  (3)),  many  assumptions  are  made  and 
several  terms  are  neglected.  Depending  on  the 
location  and  time  of  year  the  largest  of  these 
tentative  sources  of  error  might  differ  significantly. 
This  makes  finding  and  estimating  the  primary 
source  of  error  an  extremely  difficult  task  to 
perform.  Momentum  equation  error  estimates 
should  be  based  on  the  magnitudes  of  the  terms 
neglected.  However,  there  is  not  a  procedure  for 
estimating  many  of  the  neglected  terms,  thus 
requiring  estimates  based  on  the  magnitudes  of 
known  terms  and  the  assumption  that  the  un¬ 
known  terms  are  comparable.  Prior  observations 
indicate  that  the  flow  through  the  straight  is 
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Table  2 

Estimated  expected  errors  for  the  dynamics  and  measurements 
that  are  used  to  construct  the  weighing  matrix 


n 

Expected  error  ( er ) 

Continuity 

1545 

2.22  x  10-8  m/s 

Momentum 

1545 

5.50  x  10“5m2/s2 

U  &  V  smoothing 

1545 

1.65  x  10"7 1/s 

SSH  smoothing 

1545 

7.80  x  10-lol/m 

Boundary  constraints 

182 

4.87  x  10“4  m2/s 

ADCP  measurements 

12 

2.09  m2/s 

TOPEX  SSH  differences 

976 

2.18m 

n  is  the  number  of  discretized  equations  of  each  dynamical  and 
data  component  for  each  season  within  the  coefficient  matrix. 
These  values  are  required  in  Eq.  (15)  so  that  each  weighted 
component  has  an  equal  contribution  within  the  total  cost 
function  (Eq.  (12)).  Note  that  the  large  range  in  values  of 
expected  error  is  due  to  the  different  equations  and  measure¬ 
ments  having  different  physical  representation  and  thus 
different  units. 


generally  in  geostrophic  balance  (Katoh,  1993). 
Basing  error  estimates  on  either  the  Coriolis  force 
or  pressure  gradient  therefore  would  be  inap¬ 
propriate.  The  combined  magnitude  of  neglected 
terms  is  therefore  approximated  to  be  the  magni¬ 
tude  of  the  wind  stress.  Since  wind  stress  varies 
considerably  throughout  the  year,  the  average  of 
all  three  seasons  is  used  to  construct  this  expected 
error  estimate. 

The  primary  source  of  error  in  the  boundary 
constraint  is  expected  to  be  from  neglecting  the 
influence  of  river  water  outflow.  The  value  for  this 
error  is  estimated  by  considering  the  average  river 
runoff  of  South  Korea  per  unit  length  of  its 
coastline.  The  average  annual  rainfall  and  land- 
mass  area  of  South  Korea  is  1.12m/yr  and 
98,480  km2,  respectively.  If  the  percentage  of 
rainfall  that  becomes  river  runoff  is  estimated  to 
be  33.6%,  then  the  total  average  river  runoff  of 
South  Korea  is  Kr  =  1174m3/s. 

The  primary  purpose  of  the  smoothing  con¬ 
straint  is  to  filter  out  sub-grid  scale  features. 
Therefore,  the  expected  error  of  smoothing  is  set 
equal  to  the  reciprocal  of  the  square  of  the 
discretization  (1/Ax2).  Since  smoothing  is  applied 
to  both  transport  and  SSH,  the  expected  errors  are 
multiplied  by  corresponding  average  magnitude 
values  to  normalize  them.  The  average  horizontal 


transport  and  SSH  magnitude  values  are  obtained 
from  the  ADCP  and  TOPEX  data,  respectively. 

The  primary  source  of  measurement  error  is  a 
result  of  using  individual  measurements  to  repre¬ 
sent  seasonally  averaged  transports.  The  expected 
ADCP  measurement  error  is  estimated  by  comput¬ 
ing  the  RMS(t/  —  f/S)  =  8.980m2/s,  where  U 
consists  of  the  individual  transports  and  U  is 
the  seasonal-average  of  U.  Similarly,  the  repre¬ 
sentation  error  caused  by  seasonally  averaging  the 
TOPEX  data  is  RMS(?j  —  fjs)  =  9.35  cm.  Because 
the  measurements  are  seasonally  averaged  prior  to 
assimilation,  these  RMS  errors  are  divided  by  the 
equivalent  degrees  of  freedom,  \fTJx ,  where  T  is 
the  time  period  spanning  the  seasonally  averaged 
data  and  x  is  the  time  scale  of  events  (r  =  5  days; 
Teague  et  al.,  2002). 

There  are  uncertainties  in  these  expected  error 
estimates.  To  ensure  the  solution  is  not  overly 
sensitive  to  the  weights  in  the  cost  function  and 
thus  to  ensure  that  the  conclusions  reached  within 
this  study  are  well  substantiated,  each  weight  was 
individually  varied  by  plus  and  minus  one  order  of 
magnitude  to  simulate  a  wide  range  of  terms 
neglected  in  the  estimation  of  the  expected  errors. 
Comparison  of  these  solutions  with  the  solution 
presented  in  Section  5.1  reveal  only  subtle 
differences,  and  these  differences  do  not  change 
any  of  the  discussions  or  conclusions  reached  in 
this  study. 


5.  Results 

5.1.  General  circulation 

Solutions  are  obtained  for  summer,  autumn, 
and  winter  of  1999-2000  by  assimilating  the 
weighted  seasonally  averaged  ADCP  and  TOPEX 
data  with  the  weighted  barotropic  dynamics  (Fig. 
3).  The  color  shading  represents  the  SSH  solution, 
the  black  arrows  indicate  the  transport  solution, 
and  the  thick  red  arrows  indicate  the  ADCP 
measurements.  The  ADCP  measurements  are 
included  in  the  plots  as  a  reference.  It  is  important 
to  note  that  the  SSH  solutions  are  biased  since  no 
SSH  measurements  are  assimilated.  The  only  SSH 
information  that  is  included  in  the  assimilation  is 
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Fig.  3.  The  inverse  solutions  for  the  seasonal  circulation  for  (A)  summer,  (B)  autumn,  and  (C)  winter  are  obtained  by  assimilating  the 
weighted  current  meter  and  TOPEX  data  (shown  in  Fig.  2)  with  weighted  barotropic  dynamics.  The  black  vectors  represent  seasonally 
averaged  and  vertically  integrated  velocities  and  the  colored  contours  represent  SSH  (m).  ADCP  data  (red  vectors)  are  overlaid  onto 
the  figure  for  comparison. 
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the  spatial  and  seasonal  gradients  of  SSH  inferred 
from  the  pressure  gradient  term  and  TOPEX  data, 
respectively. 

In  all  three  seasons,  the  solution  shows  inflow 
through  the  open  boundaries  of  the  domain 
southwest  of  the  strait  with  some  flow  entering 
through  the  Cheju  Strait.  The  incoming  water 
converges  to  the  center  of  the  strait  just  south  of 
Tsushima  Island.  The  flow  then  branches  through 
the  channels  east  and  west  of  Tsushima  Island,  and 
upon  exiting  the  strait,  forms  the  EKWC,  the 
NSB,  and  the  OSB.  There  are  significant  differ¬ 
ences  in  the  branching  of  the  strait  outflow 
between  the  three  seasonal  solutions.  The  summer 
and  autumn  transports  through  the  eastern  chan¬ 
nel  are  similar,  whereas  the  transport  through  the 
western  channel  increases  substantially  from  sum¬ 
mer  to  autumn.  During  winter,  the  transport 
through  both  channels  decreases  considerably.  In 
the  autumn  solution,  there  is  a  small  cyclonic  eddy 
just  northeast  of  Tsushima  Island.  In  the  winter 
solution,  there  is  a  small  anticyclonic  eddy  south¬ 
west  of  the  eastern  channel.  This  eddy  is  a  result  of 
there  being  no  data  at  station  S5  and  the  velocity 
at  station  S6  being  directed  in  the  southwestward 
direction. 

The  NSB  and  OSB  are  very  distinct  in  the 
solutions  for  all  three  seasons.  At  the  exit  of  the 
western  channel  (about  35.5°  latitude),  a  majority 
of  the  transport  branches  toward  the  southeast 
following  the  southern  ridge  of  the  Ullung  basin 
and  contributes  to  the  transport  within  the  OSB. 
This  diversion  of  water  to  the  OSB  causes  the 
EKWC  to  be  much  weaker  than  has  been 
observed.  The  dynamical  equation  residuals  re¬ 
sulting  from  these  solutions  provide  ancillary 
information  needed  in  order  to  understand  this 
result,  which  is  discussed  in  Section  6. 

5.2.  Convergence  of  solution 

In  order  to  determine  the  level  of  dependence 
the  solution  has  on  just  the  ADCP  measurements 
and  how  well  these  measurements  constrain  the 
system,  the  assimilation  of  TOPEX  data  is 
removed  and  the  ADCP  measurements  from  the 
two  lines  of  moorings  are  separately  removed. 
Comparison  of  autumn  solutions  using  only  the 


northern  and  then  only  the  southern  lines  of 
moorings  (Figs.  4B  and  C,  respectively)  with  the 
autumn  solution  that  uses  all  moorings  (Fig.  4A) 
reveals  that  even  with  fewer  measurements  the 
solution  remains  fairly  consistent.  The  spatially 
averaged  percent  difference  between  transport 
solutions  that  use  just  the  northern  line  of 
moorings  and  all  moorings  is  roughly  20%. 
Likewise,  the  difference  using  just  the  southern 
line  of  moorings  is  roughly  30%.  Therefore,  the 
solution  using  just  the  northern  line  of  measure¬ 
ments  is  able  to  reproduce  the  full  solution  more 
accurately  than  using  just  the  southern  line.  The 
reason  for  this  result  is  that  there  is  a  large 
discrepancy  between  the  dynamics  and  measure¬ 
ments  at  the  exit  of  the  western  channel,  which  is 
in  the  northern  half  of  the  domain.  The  cause  for 
this  discrepancy  will  be  analyzed  in  Section  5.3. 

Comparison  of  the  autumn  solutions  with  and 
without  the  assimilation  of  TOPEX  data  (Figs.  3B 
and  4A,  respectively)  reveal  that  the  basic  struc¬ 
ture  of  the  flow  fields  is  roughly  the  same.  There  is, 
however,  a  significant  change  in  the  flow  field 
around  Cheju  Island.  The  addition  of  TOPEX 
data  causes  a  large  reduction  in  the  CWC  and  a 
large  increase  in  northward  flow  west  of  Cheju 
Island.  Also,  this  comparison  shows  a  significant 
difference  in  the  magnitudes  of  SSH.  This  is  to  be 
expected  since  the  SSH  solution  without  the 
assimilation  of  TOPEX  data  is  only  a  function  of 
the  pressure  gradient  term.  By  including  measure¬ 
ments  of  seasonal  SSH  difference  in  the  assimila¬ 
tion,  the  SSH  is  shifted  in  order  for  all  seasonal 
solutions  to  satisfy  these  data.  The  inclusion  of 
TOPEX  data  also  provides  additional  spatial 
coverage  that  helps  constrain  the  solution  away 
from  the  ADCP  mooring  locations.  Transport 
solutions  away  from  the  ADCP  measurements 
have  a  higher  magnitude  and  smaller  spatial  scales, 
which  correspond  to  steeper  SSH  gradients. 

Even  though  the  TOPEX  data  provide  useful 
information  to  the  assimilation  system,  it  is  not  the 
primary  reason  they  were  included.  The  main 
reason  for  this  inclusion  is  that  they  help  constrain 
the  weighted  least  squares  problem  and  they 
improve  the  condition  number  of  the  normal 
matrix.  To  demonstrate  that  adding  TOPEX  data 
improves  the  condition  number,  the  values  of  the 
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Fig.  4.  Inverse  solutions  for  autumn  obtained  by  assimilating  just  the  weighted  current  meter  data  (same  as  Fig.  3b  but  without  the 
assimilation  of  TOPEX  data).  Comparison  of  inverse  solutions  resulting  from  the  assimilation  of  (A)  all,  (B)  the  northern  line,  and  (C) 
the  southern  line  of  current  meters  shows  that  observations  match  the  barotropic  dynamics  reasonably  well  at  the  locations  of  the 
missing  current  meters. 
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Fig.  5.  The  minimization  of  the  seasonal  components  of  the 
total  cost  function  (including  the  assimilation  of  TOPEX  data) 
as  it  iterates  through  the  conjugate  gradient.  The  displayed 
autumn  solution  without  the  assimilation  of  TOPEX  data 
shows  that  reducing  the  number  of  observations  increases  the 
number  of  iterations  required  for  the  cost  function  to  converge. 


cost  function  (Eq.  (12))  after  each  iteration  of  the 
conjugate  gradient  are  displayed  for  each  season 
(Fig.  5).  The  convergence  of  the  cost  function  for 
the  autumn  solution  without  the  assimilation  of 
TOPEX  data  is  also  displayed  for  comparison. 
These  curves  provide  a  measure  of  how  well  the 
state  satisfies  all  of  the  contributors  to  the  total 
cost  function,  which  include  errors  to  the  dy¬ 
namics,  constraints,  and  measurements.  As  dis¬ 
cussed  in  Section  4.3,  the  weights  for  each 
contributor  are  divided  by  its  number  of  equations 
(the  number  of  grid  points,  boundary  points,  or 
measurements)  such  that  each  contributor  is 
considered  equally  during  the  convergence  of  the 
cost  function.  This  means  that  if  each  discretized 
equation  from  all  of  the  contributors  is  exactly 
satisfied  within  expected  error,  then  the  total  cost 
would  be  equal  to  the  number  of  contributors 
(nine).  The  cost  functions  in  Fig.  5  indicate  that 
the  converged  solution  for  each  season  satisfies  the 
overall  contribution  from  dynamics,  constraints, 
and  measurements  to  within  expected  error.  The 
total  residual  for  the  summer  and  winter  solutions 
is  within  10%  of  the  expected  error.  This  plot  also 
shows  that  the  solution  for  winter  satisfies  the 
overall  specified  barotropic  dynamics  the  best,  and 
the  autumn  solution  contains  the  largest  discre¬ 
pancies  between  the  dynamics  and  data. 

However,  even  though  the  solutions  are  con¬ 
verged,  the  number  of  iterations  that  are  required 


to  achieve  this  convergence  is  more  than  optimal. 
As  mentioned  in  Section  4.2,  one  to  three  loops 
through  the  conjugate  gradient  routine  is  typically 
sufficient  to  obtain  a  converged  solution  for  a  well- 
conditioned  system,  where  each  loop  consists  of 
the  number  of  iterations  equal  to  the  size  of  the 
state.  The  autumn  solution  without  the  assimila¬ 
tion  of  TOPEX  data  requires  250  loops  (roughly 
1.1  million  iterations)  through  the  conjugate 
gradient  to  converge  (Fig.  5).  Even  after  this  many 
iterations,  the  total  cost  still  has  a  slight  downward 
trend,  indicating  the  problem  is  poorly  condi¬ 
tioned  without  the  assimilation  of  TOPEX  data. 
Whereas,  the  autumn  solution  with  the  assimila¬ 
tion  of  TOPEX  data  only  requires  75  loops 
(roughly  1  million  iterations)  through  the  con¬ 
jugate  gradient  in  order  to  converge.  Note  that  the 
inclusion  of  TOPEX  data  requires  all  three 
seasonal  solutions  to  be  solved  simultaneously. 
Therefore,  the  size  of  the  state  and  hence  the 
number  of  iterations  contained  in  one  loop  is  three 
times  larger  than  that  without  the  assimilation  of 
TOPEX  data.  In  other  words,  the  total  number  of 
iterations  required  to  solve  all  three  seasonal 
solutions  with  the  assimilation  of  TOPEX  data  is 
less  than  that  required  to  solve  one  season  without 
TOPEX  data. 


5.3.  Dynamic  residuals 

Since  the  computed  solutions  are  a  best  fit  to  the 
equations  of  motion,  constraints,  and  measure- 

Table  3 

Spatially  averaged  residuals  of  each  component  of  the 
barotropic  dynamics  and  data  for  each  season  divided  by 
corresponding  expected  errors 


Summer 

Autumn 

Winter 

Continuity 

5.317e-4 

8.761e-4 

5.456e-4 

Momentum  ' 

0.5256 

0.9107 

0.5791 

U  &  V  smoothing 

0.2174 

0.3076 

0.2074 

SSH  smoothing 

0.0526 

0.0771 

0.0517 

Boundary  constraints 

1.620e-4 

2.321e-4 

1.508e-4 

ADCP  measurements 

0.0420 

0.0539 

0.0403 

TOPEX  SSH  differences 

0.6650 

0.9525 

The  units  of  these  residuals  are  in  number  of  standard 
deviations  (Std)  of  the  corresponding  expected  errors. 
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Fig.  6.  The  momentum  equation  residuals  show  where  during  the  (A)  summer,  (B)  autumn,  and  (C)  winter  seasons  that  the  assumed 
barotropic  dynamics  are  in  error.  These  residuals  are  in  units  of  number  of  standard  deviations  (Std)  relative  to  the  expected  error  and 
are  considered  significant  if  they  are  greater  than  3  Std.  The  largest  residuals  occur  in  autumn  at  the  exit  of  the  western  channel.  It  is  in 
this  area  that  baroclinic  dynamics  are  expected  to  be  important  based  on  previous  studies. 
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ments,  the  solution  is  going  to  satisfy  none  of  these 
exactly  (Eq.  (11)).  The  dynamic  and  measurement 
residuals  are  computed  and  examined  in  order  to 
determine  how  well  the  solutions  satisfy  the 
assumed  barotropic  dynamics.  The  absolute  value 
of  these  residuals  is  divided  by  the  corresponding 
expected  errors  (|Ax  —  h\j (er))  to  identify  areas  in 
which  residuals  are  much  larger  than  expectations. 
These  resulting  measures  of  error  are  therefore  in 
terms  of  number  of  standard  deviations  (Std)  of 
the  expected  error.  Spatial  averages  of  the  normal¬ 
ized  residuals  are  displayed  in  Table  3  for  each 
component  and  season,  and  give  an  overall 
indication  of  how  well  the  solution  satisfies  the 
prescribed  dynamics  and  measurements.  In  gen¬ 
eral,  the  spatially  averaged  normalized  residuals 
are  roughly  equivalent  for  summer  and  winter, 
whereas  for  autumn  these  residuals  are  higher.  The 
mean  residuals  show  that  continuity  and  the 
boundary  constraints  are  satisfied  well  within 
expected  values,  residuals  in  the  smoothing  of 
SSH  and  the  ADCP  data  are  within  0.1  Std, 
residuals  in  the  smoothing  of  horizontal  transport 
are  within  0.5  Std,  and  agreement  to  the  momen¬ 
tum  equations  and  the  TOPEX  measurements  are 
within  1  Std.  The  distribution  of  the  normalized 
barotropic  momentum  equation  residuals  (Fig.  6) 
indicates  that  the  solution  satisfies  the  assumed 
dynamics  to  within  1  Std  of  expected  values  over 
most  of  the  domain,  especially  during  winter. 
There  are,  however,  small  pockets  of  relatively 
large  error  at  some  of  the  ADCP  mooring 
locations  during  all  three  seasons,  and  a  large 
area  of  error  at  the  exit  of  the  western  channel 
during  autumn.  At  these  locations,  momentum 
equation  residuals  exceed  3  Std  of  the  expected 
error. 

It  is  believed  that  the  significantly  large  mo¬ 
mentum  residual  at  the  exit  of  the  western 
channel  during  autumn  is  a  result  of  the  barotropic 
momentum  equations  not  being  able  to  resolve 
baroclinic  processes.  Observations  have  been 
made  in  previous  studies  that  suggest  that 
strong  baroclinic  processes  occur  in  this  region 
and  peak  during  this  season  (Isobe,  1995).  To  test 
this  theory,  solutions  are  computed  by  assimilating 
the  ADCP  and  TOPEX  measurements  with  the 
2.5-layered  baroclinic  system  of  equations  dis- 


Table  4 

Spatially  averaged  residuals  of  each  component  of  the 
baroclinic  dynamics  and  data  for  each  season  and  layer  divided 
by  corresponding  expected  errors 


Summer 

Autumn 

Winter 

Continuity 

Layer  1 

5.856e-4 

4.130e-4 

4.121e-4 

Layer  2 

0.0972 

0.1373 

0.0916 

Momentum 

Layer  1 

0.3538 

0.5346 

0.3968 

Layer  2 

0.4620 

0.5665 

0.3939 

U  &  V  smoothing 

Layer  1 

0.2096 

0.2847 

0.1805 

Layer  2 

0.3304 

0.4175 

0.2389 

r\  smoothing 

Layer  1 

0.0674 

0.0860 

0.0532 

Layer  2 

0.2769 

0.3559 

0.2221 

Layer  3 

0.0504 

0.0542 

0.0362 

Boundary  constraints 

Layer  1 

1.240e-4 

L591e-4 

1.360e-4 

Layer  2 

1.577e-4 

1.969e-4 

1.752e-4 

ADCP  measurements 

Layer  1 

0.0220 

0.0486 

0.0276 

Layer  2 

0.0307 

0.0283 

0.0182 

TOPEX  SSH  differences 

0.3702 

0.5002 

The  units  of  these  errors  are  in  number  of  standard  deviations 
(Std)  of  the  corresponding  expected  errors. 


cussed  in  Section  3.2.  The  spatially  averaged 
residuals  relative  to  their  expected  values 
resulting  form  this  assimilation  are  presented  in 
Table  4  for  each  component,  season,  and  layer. 
Comparison  of  these  residuals  with  those  resulting 
from  using  the  barotropic  system  of  equations 
(Table  3)  reveals  that  the  mean  residuals  of 
continuity,  smoothing,  and  the  boundary  con¬ 
straint  have  comparable  values.  The  mean 
residuals  of  the  momentum  equations,  ADCP 
data,  and  TOPEX  data,  however,  are  roughly  half 
in  the  2.5-layered  baroclinic  system  compared  to 
the  corresponding  residuals  in  the  barotropic 
system.  Thus,  the  baroclinic  momentum  equations 
are  in  better  agreement  with  the  measurements 
than  the  barotropic  momentum  equations.  By 
comparing  the  distribution  of  momentum  resi¬ 
duals  during  autumn  between  these  two  systems 
of  equations  (Fig.  7),  it  is  apparent  that  the 
inclusion  of  baroclinic  processes  in  the  momentum 
equations  reduces  the  residual  at  the  exit  of  the 
western  channel. 
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6.  Discussion  and  conclusions 

Velocity  and  SSH  measurements  are  assimilated 
with  both  a  system  of  single-layered  barotropic 
and  2.5-layered  baroclinic  equations  for  the 
summer,  autumn,  and  winter  seasons  of 
1999-2000  using  a  variational  assimilation  ap¬ 
proach.  Since  the  weighted  data  and  dynamics 
cannot  be  expected  to  match  each  other  exactly, 
the  optimal  solution  will  cause  residuals  in  both. 
The  barotropic  momentum  equation  residuals 
(Fig.  6)  reveal  that  the  best-fit  solution  does  not 
satisfy  the  assumed  dynamics  within  expected 
values  at  the  exit  of  the  western  channel  during 
autumn.  The  location  and  season  of  this  large 
momentum  residual  coincide  with  the  region  and 
time  of  year  that  have  been  observed  in  previous 
studies  to  exhibit  baroclinic  dynamics.  For  exam¬ 
ple,  Isobe  (1995)  assimilates  ADCP-CTD  survey 
data  into  a  5-layered  model  for  the  summer  and  a 
2-layered  model  for  the  winter,  and  demonstrates 
that  the  observed  seasonal  variability  in  the  sea- 
level  difference  across  the  Tsushima/Korean  Strait 
just  north  of  Tsushima  Island  largely  contains  the 
baroclinic  motion  caused  by  the  Bottom  Cold 
Water  (BCW),  which  is  believed  to  originate  in  the 
Japan/East  Sea  and  penetrate  almost  exclusively 
into  the  western  channel.  Isobe  (1995)  also 
concludes  that  the  sea  level  differences  across  the 
strait  is  8  cm  in  the  summer  and  3  cm  in  the  winter, 
which  corresponds  with  the  amount  of  BCW 
observed  in  the  western  channel  being  at  its 
maximum  during  summer  and  minimum  during 
winter. 

In  a  more  recent  hydrographic  study,  Cho  and 
Kim  (2000)  observed  the  BCW  moving  up  onto 
the  shelf  of  the  western  channel  during  summer 
and  retreating  back  to  the  Ullung  basin  during 
winter.  In  addition  to  showing  the  nonexistence  of 
the  BCW  during  winter,  observations  were  made 
that  indicate  that  the  EKWC  is  also  absent  during 
this  season.  Cho  and  Kim  (2000)  explain  that 
when  the  lower  layer  is  absent,  the  northward  flow 
in  the  upper  layer  has  a  positive  relative  vorticity 
due  to  the  increasing  depth  from  south  to  north, 
therefore  causing  the  flow  exiting  the  western 
channel  to  veer  eastward.  In  contrast  during 
summer,  the  existence  of  the  BCW  causes  the 


upper  layer  to  shrink  as  the  flow  moves  northward, 
therefore  creating  negative  relative  vorticity  and 
causing  the  EKWC  to  flow  along  the  western 
boundary. 

In  the  barotropic  results  (Fig.  3),  the  EKWC  is 
virtually  nonexistent  during  summer  and  winter 
and  weaker  than  observed  during  autumn.  Since 
the  barotropic  dynamics  do  not  account  for  the 
baroclinic  effects  of  the  BCW,  then  according  to  Cho 
and  Kim’s  (2000)  hypothesis  it  is  consistent  that  the 
flow  exiting  the  western  channel  in  these  results  veers 
eastward  instead  of  forming  the  EKWC.  The 
assimilation  of  measurements  with  a  2.5-layered 
baroclinic  system  of  equations  shows  that  the 
incorporation  of  baroclinic  effects  results  in  a 
significant  reduction  in  momentum  error  at  the  exit 
of  the  western  channel  during  autumn  (Fig.  7).  This 
implies  that  there  is  a  failure  of  the  barotropic 
approximation  within  this  region.  Even  though  it  is 
widely  believed  that  the  flow  through  the  Tsushima 
Strait  is  mostly  barotropic,  the  results  from  this 
study  indicate  that  some  other  processes  are 
prevalent  at  the  exit  of  the  western  channel  and  are 
capable  of  seriously  affecting  the  overall  dynamics  of 
the  flow.  The  most  plausible  source  of  these 
processes  within  this  area  is  baroclinic  forcing. 
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